function P = ellipsefit(X,Y)
x=X;
y=Y;
xy=x.*y;
x2=x.*x;
y2=y.*y;
x3=x2.*x;
y3=y2.*y;
x2y=x2.*y;
xy2=x.*y2;
x4=x2.*x2;
y4=y2.*y2;
x3y=x3.*y;
xy3=x.*y3;
x2y2=x2.*y2;
M=[ sum(x4)    sum(x3y)    sum(x2y2)    sum(x3)    sum(x2y);
    sum(x3y)   sum(x2y2)   sum(xy3)     sum(x2y)   sum(xy2);
    sum(x2y2)  sum(xy3)    sum(y4)      sum(xy2)   sum(y3) ;
    sum(x3)    sum(x2y)    sum(xy2)     sum(x2)    sum(xy) ;
    sum(x2y)   sum(xy2)    sum(y3)      sum(xy)    sum(y2) ];
N = [sum(x2) sum(xy) sum(y2) sum(x) sum(y)]';
P=M\N;